library(tidyverse)
library(plotrix) #std.error
library(hypr) #for generating contrast coding
library(ggthemes)
library(ggh4x)
library(dagitty) #dagitty for causal diagrams
library(DiagrammeR) #grViz for causal diagrams
library(DiagrammeRsvg)
library(rsvg)
library(plotly) #ggplotly for interactive plots
library(brms) #bayesian regression
library(bayestestR) #p_direction, hdi
library(tidybayes) #add_epred_draws
library(emmeans)
library(doParallel) #set the number of cores
library(ppcor) #partial correlation
# # for power analysis
# library(lme4)
# library(mixedpower)
### Set global options
options(digits = 3) # set the default number of digits to 3
cores = as.integer(detectCores(logical = FALSE) * 0.7) # set the number of cores to use
registerDoParallel(cores=cores) # register the number of cores to use for parallel processing
options(mc.cores = cores)
options(brms.backend = "cmdstanr") #this will speed up the model fitting
### MCMC options
niter = 20000 #number of iterations
nwu = 2000 #number of warmups
### Rmd settings
knitr::opts_chunk$set(echo=TRUE, warning=FALSE, message=FALSE)
### Custom functions
########### pp_check_each_round ############
pp_check_each_round <- function(m, data, round_info) {
df_sub = filter(data, round == round_info)
p = pp_check(m,
type = "bars",
ndraws = 100,
newdata = df_sub) +
coord_cartesian(xlim = c(0, 10)) +
ggtitle(round_info)
return(p)
}
########### pp_check_each_condition ############
pp_check_each_condition <- function(m, data, condition_info) {
df_sub = filter(data, condition == condition_info)
p = pp_check(m,
type = "bars",
ndraws = 100,
newdata = df_sub) +
coord_cartesian(xlim = c(0, 10)) +
ggtitle(condition_info)
return(p)
}
########### plot_posterior ############
plot_posterior <- function(model, model2=NA, model3=NA,
interaction=FALSE, include_intercept=FALSE,
xlim_cond=1.5, xlim_round=2, xlim_lex=0.15){
### extract the posterior draws
posterior_beta1 <- model %>%
gather_draws(`b_.*`, regex = TRUE) %>%
mutate(intercept = str_detect(.variable, "Intercept"),
component = ifelse(str_detect(.variable, ":"), "Interaction",
ifelse(str_detect(.variable, "round"), "Round",
ifelse(str_detect(.variable, "Intercept"), "Intercept",
ifelse(str_detect(.variable, "lex_align"), "N. lex align",
"Visibility")))))
if (length(model2) == 1){ #if model2 is NA
posterior_beta = posterior_beta1
} else {
posterior_beta1 <- posterior_beta1 %>%
filter(.variable != "b_num_lex_align")
posterior_beta2 <- model2 %>%
gather_draws(`b_.*`, regex = TRUE) %>%
filter(.variable == "b_num_lex_align") %>%
mutate(component = "N. lex align")
posterior_beta3 <- model3 %>%
gather_draws(`b_.*`, regex = TRUE) %>%
filter(.variable == "b_condition_sumAO_Sym") %>%
mutate(component = "Visibility")
posterior_beta = rbind(posterior_beta1, posterior_beta2, posterior_beta3)
}
posterior_beta = posterior_beta %>%
mutate(.variable = recode(.variable,
"b_Intercept" = "Intercept",
"b_conditionAsymAV" = "SymAV--AsymAV",
"b_conditionAO" = "SymAV--AO",
"b_conditionAsym_Sym" = "AsymAV--SymAV",
"b_conditionAO_Asym" = "AO--AsymAV",
"b_condition_sumAO_Sym" = "AO--SymAV",
"b_num_lex_align" = "N. lex align",
"b_round_c" = "Centered round",
"b_log_round_c" = "Centered log(round)",
"b_roundR12" = "R1--R2",
"b_roundR23" = "R2--R3",
"b_roundR34" = "R3--R4",
"b_roundR45" = "R4--R5",
"b_roundR56" = "R5--R6",
"b_conditionAsymAV:round1" = "SymAV--AsymAV: R1--R2",
"b_conditionAsymAV:round2" = "SymAV--AsymAV: R2--R3",
"b_conditionAsymAV:round3" = "SymAV--AsymAV: R3--R4",
"b_conditionAsymAV:round4" = "SymAV--AsymAV: R4--R5",
"b_conditionAsymAV:round5" = "SymAV--AsymAV: R5--R6",
"b_conditionAO:round1" = "SymAV--AO: R1--R2",
"b_conditionAO:round2" = "SymAV--AO: R2--R3",
"b_conditionAO:round3" = "SymAV--AO: R3--R4",
"b_conditionAO:round4" = "SymAV--AO: R4--R5",
"b_conditionAO:round5" = "SymAV--AO: R5--R6",
"b_conditionAsym_Sym:round_c" = "Centered round: Asym--Sym",
"b_conditionAO_Sym:round_c" = "Centered round: AO--Sym",
"b_conditionAO_Asym:round_c" = "Centered round: AO--Asym",
"b_conditionAsym_Sym:log_round_c" = "Centered log(round): Asym--Sym",
"b_conditionAO_Sym:log_round_c" = "Centered log(round): AO--Sym",
"b_conditionAO_Asym:log_round_c" = "Centered log(round): AO--Asym"),
.variable = factor(.variable,
levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV",
"R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
"N. lex align")),
component = factor(component,
levels = c("Intercept", "Visibility", "Round",
"Interaction", "N. lex align")))
if (include_intercept == F){
posterior_beta = posterior_beta %>% filter(component != "Intercept")
}
### change variables if only main effects are plotted
if (interaction == F) {
posterior_beta = filter(posterior_beta, !str_detect(.variable, ":"))
fill_manual_values = c("steelblue", "steelblue", "steelblue")
} else{
fill_manual_values = c("steelblue", "steelblue", "steelblue", "steelblue")
}
### plot the posterior distributions
p_posterior = ggplot(posterior_beta,
aes(x = .value, y = fct_rev(.variable),
fill = component)) +
geom_vline(xintercept = 0, size = 1) +
stat_halfeye(aes(slab_alpha = intercept),
normalize = "panels",
slab_alpha = .5,
.width = c(0.89, 0.95),
point_interval = "median_hdi") +
scale_fill_manual(values = fill_manual_values) +
scale_slab_alpha_discrete(range = c(0.8, 0.4)) +
guides(fill = "none", slab_alpha = "none") +
labs(x = "Coefficient", y = "Effect") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
strip.background = element_blank(),
panel.grid.major.x = element_line(color = "grey90",
linetype = "solid",
size = 0.5),
panel.grid.major.y = element_line(color = "grey90",
linetype = "solid",
size = 0.5),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
facet_wrap(vars(component), ncol = 3, scales = "free") +
facetted_pos_scales(
x = list(
scale_x_continuous(limits = c(-xlim_cond, xlim_cond)),
scale_x_continuous(limits = c(-xlim_round, xlim_round)),
scale_x_continuous(limits = c(-xlim_lex, xlim_lex),
breaks = c(-0.1, 0, 0.1))))
return(p_posterior)
}
########### pp_update_plot ############
pp_update_plot <- function(post_sample, model_type="nb", interaction=TRUE){
sum = ifelse("b_condition_sumAO_Sym" %in% colnames(post_sample), T, F)
round_bd = ifelse("b_roundR12" %in% colnames(post_sample), T, F)
intercept = ggplot(post_sample) +
geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Intercept') +
theme_classic()
### Visibility condition
if (sum == F){
cond1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
cond2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO--Asym') +
theme_classic()
} else {
cond1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_condition_sumAO_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO--Sym') +
theme_classic()
cond2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_condition_sumAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
}
### Round
if (interaction) {
if (round_bd){
r1 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR12), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R1--R2') +
theme_classic()
r2 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR23), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R2--R3') +
theme_classic()
r3 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR34), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R3--R4') +
theme_classic()
r4 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR45), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R4--R5') +
theme_classic()
r5 = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_roundR56), fill="#FC4E07", color="black",alpha=0.6) +
xlab('R5--R6') +
theme_classic()
} else {
if ("b_round_c" %in% colnames(post_sample)) {
round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered round') +
theme_classic()
if (sum == F){
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: Asym--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAO_Asym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: AO--Asym') +
theme_classic()
} else {
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAO_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: AO--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAsym_Sym:round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered Round: Asym--Sym') +
theme_classic()
}
} else {
round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round)') +
theme_classic()
if (sum == F){
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): Asym--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_conditionAO_Asym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): AO--Asym') +
theme_classic()
} else {
cond1_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAO_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): AO--Sym') +
theme_classic()
cond2_round = ggplot(post_sample) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(`b_condition_sumAsym_Sym:log_round_c`), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round): Asym--Sym') +
theme_classic()
}
}
}
}
if (model_type == "nb"){
shape = ggplot(post_sample) +
geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Shape') +
scale_x_continuous(limits = c(0, 10)) +
theme_classic()}
else if (model_type == "zinb") {
shape = ggplot(post_sample) +
geom_density(aes(prior_shape), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(shape), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Shape') +
scale_x_continuous(limits = c(0, 10)) +
theme_classic()
zi = ggplot(post_sample) +
geom_density(aes(prior_zi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(zi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Zero-inflation') +
theme_classic()}
else if (model_type == "zibt"){
phi = ggplot(post_sample) +
geom_density(aes(prior_phi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(phi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Precision') +
theme_classic()
zoi = ggplot(post_sample) +
geom_density(aes(prior_zoi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(zoi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Zero-inflation') +
theme_classic()
coi = ggplot(post_sample) +
geom_density(aes(prior_coi), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(coi), fill="#FC4E07", color="black",alpha=0.6) +
xlab('One-inflation') +
theme_classic()
}
### display the plots
if (interaction==F){
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=2)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
}
} else {
if (round_bd == T){
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, ncol=3)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, phi, zoi, coi, ncol=3)
}
} else {
if (model_type=="nb"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, shape, ncol=3)}
else if (model_type=="zinb"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, shape, zi, ncol=3)}
else if (model_type=="zibt"){
gridExtra::grid.arrange(intercept, cond1, cond2, round,
cond1_round, cond2_round, phi, zoi, coi, ncol=3)
}
}
}
}
### trial_info.csv
df_trial_info = read.csv("data/trial_info.csv", stringsAsFactors = TRUE) %>%
filter(num_words != 0) %>% # remove trials that were accidentally skipped) %>%
mutate(pair = factor(pair),
round_c = as.integer(round) - mean(as.integer(round)),
round_n = factor(round),
round = factor(paste0("R", round)),
log_round = log(as.integer(round)),
log_round_c = log_round - mean(log_round),
condition = factor(condition,
levels = c("Sym", "Asym", "AO"),
labels = c("SymAV", "AsymAV", "AO")),
condition_sum = condition,
duration_s = duration_ms / 1000,
n_words_100 = num_words / 100,
n_iconic_per_100words = num_iconic_gestures / n_words_100,
n_iconic_100 = num_iconic_gestures / 100) %>%
dplyr::select(pair, condition, round, round_n,
round_c, log_round, log_round_c,
duration_ms, duration_s, duration_m,
director:n_iconic_100) %>%
rename(trial_duration_s = duration_s, trial_duration_m = duration_m)
head(df_trial_info)
### condition info
df_condition_info = read.csv("data/condition_info.csv", stringsAsFactors = TRUE) %>%
mutate(pair = factor(pair),
condition = factor(condition,
levels = c("Sym", "Asym", "AO"),
labels = c("SymAV", "AsymAV", "AO")))
### participants
df_participant = read_tsv("data/participant_info.txt") %>%
mutate(pair = factor(as.integer(pair)),
condition = factor(condition))
df_participant_used = df_participant %>%
filter(used == "Y") %>%
mutate(condition = factor(condition, levels = c("Sym", "Asym", "AO")))
df_participant_used
mean_age = mean(df_participant_used$age)
sd_age = sd(df_participant_used$age)
range_age = paste(range(df_participant_used$age)[1], range(df_participant_used$age)[2], sep = " - ")
n_participants = nrow(df_participant_used)
n_female = nrow(filter(df_participant_used, gender == "F"))
n_male = nrow(filter(df_participant_used, gender == "M"))
data.frame(mean_age, sd_age, range_age, n_participants, n_female, n_male)
### dyads
# make a dataframe with dyad information (for each pair, gender & age of speaker 1 & 2)
df_dyad = df_participant_used %>%
dplyr::select(-notes) %>%
pivot_wider(names_from = speaker,
values_from = c("gender", "age")) %>%
mutate(mixed_gender = ifelse(gender_A == gender_B, "N", "Y"),
mixed_gender = ifelse(is.na(mixed_gender), "Y", mixed_gender),
gender_dyad = paste0(gender_A, "_", gender_B))
df_dyad %>%
group_by(condition, mixed_gender) %>%
summarize(n = n())
df_dyad %>%
group_by(condition, mixed_gender, gender_dyad) %>%
summarize(n = n())
summarize_demographics <- function(df) {
df %>%
summarize(total_s = sum(trial_duration_s),
total_m = total_s / 60,
### trial duration ###
mean_trial_dur_s = mean(trial_duration_s),
mean_trial_dur_m = mean(trial_duration_m),
sd_trial_dur_m = sd(trial_duration_m),
lci_trial_dur_m = mean_trial_dur_m - 1.96 * sd_trial_dur_m / sqrt(n()),
uci_trial_dur_m = mean_trial_dur_m + 1.96 * sd_trial_dur_m / sqrt(n()),
### words ###
# number of words
num_words_total = sum(num_words),
mean_num_words = mean(num_words),
num_words_100 = mean(num_words / 100),
sd_num_words = sd(num_words),
lci_num_words = mean_num_words - 1.96 * sd_num_words / sqrt(n()),
uci_num_words = mean_num_words + 1.96 * sd_num_words / sqrt(n()),
num_words_per_min = num_words_total / total_m,
# number of content words
num_content_total = sum(num_content_words),
mean_num_content = mean(num_content_words),
sd_num_content = sd(num_content_words),
lci_num_content = mean_num_content - 1.96 * sd_num_content / sqrt(n()),
uci_num_content = mean_num_content + 1.96 * sd_num_content / sqrt(n()),
num_content_per_min = num_content_total / total_m,
### lexical alignment ###
# raw frequency
num_lex_align_total = sum(num_lex_align, na.rm = T),
mean_num_lex_align = mean(num_lex_align, na.rm = T),
sd_num_lex_align = sd(num_lex_align, na.rm = T),
lci_num_lex_align = mean_num_lex_align - 1.96 * sd_num_lex_align / sqrt(n()),
uci_num_lex_align = mean_num_lex_align + 1.96 * sd_num_lex_align / sqrt(n()),
num_lex_align_per_min = num_lex_align_total / total_m,
# rate per 100 words
mean_lex_align_per_100words = mean(num_lex_align / num_words_100, na.rm=T),
sd_lex_align_per_100words = sd(num_lex_align / num_words_100, na.rm=T),
se_lex_align_per_100words = std.error(num_lex_align / num_words_100, na.rm=T),
lci_lex_align_per_100words = mean_lex_align_per_100words - 1.96 * se_lex_align_per_100words,
uci_lex_align_per_100words = mean_lex_align_per_100words + 1.96 * se_lex_align_per_100words,
### iconic gestures ###
# raw frequency
num_iconic_total = sum(num_iconic_gestures, na.rm = T),
mean_num_iconic = mean(num_iconic_gestures, na.rm = T),
sd_num_iconic = sd(num_iconic_gestures, na.rm = T),
lci_num_iconic = mean_num_iconic - 1.96 * sd_num_iconic / sqrt(n()),
uci_num_iconic = mean_num_iconic + 1.96 * sd_num_iconic / sqrt(n()),
num_iconic_per_min = num_iconic_total / total_m,
# rate per 100 words
mean_iconic_per_100words = mean(num_iconic_gestures / num_words_100, na.rm=T),
sd_iconic_per_100words = sd(num_iconic_gestures / num_words_100, na.rm=T),
se_iconic_per_100words = std.error(num_iconic_gestures / num_words_100, na.rm=T),
lci_iconic_per_100words = mean_iconic_per_100words - 1.96 * se_iconic_per_100words,
uci_iconic_per_100words = mean_iconic_per_100words + 1.96 * se_iconic_per_100words,
### gestural alignment ###
# raw frequency
num_gest_align_total = sum(num_gestural_align, na.rm = T),
mean_num_gest_align = mean(num_gestural_align, na.rm = T),
sd_num_gest_align = sd(num_gestural_align, na.rm = T),
lci_num_gest_align = mean_num_gest_align - 1.96 * sd_num_gest_align / sqrt(n()),
uci_num_gest_align = mean_num_gest_align + 1.96 * sd_num_gest_align / sqrt(n()),
num_gest_align_per_min = num_gest_align_total / total_m,
# rate per 100 words
mean_gest_align_per_100words = mean(num_gestural_align / num_words_100, na.rm=T),
sd_gest_align_per_100words = sd(num_gestural_align / num_words_100, na.rm=T),
se_gest_align_per_100words = std.error(num_gestural_align / num_words_100, na.rm=T),
lci_gest_align_per_100words = mean_gest_align_per_100words - 1.96 * se_gest_align_per_100words,
uci_gest_align_per_100words = mean_gest_align_per_100words + 1.96 * se_gest_align_per_100words,
# rate per iconic gestures
mean_gest_align_prop = mean(num_gestural_align / num_iconic_gestures, na.rm=T),
sd_gest_align_per_prop = sd(num_gestural_align / num_iconic_gestures, na.rm=T),
se_gest_align_per_prop = std.error(num_gestural_align / num_iconic_gestures, na.rm=T),
lci_gest_align_per_prop = mean_gest_align_prop - 1.96 * se_gest_align_per_prop,
uci_gest_align_per_prop = mean_gest_align_prop + 1.96 * se_gest_align_per_prop,
### number of trials ###
trial_n = n()) %>%
ungroup()
}
summarize_dur <- function(df){
df %>%
summarize(mean_dur_m = mean(total_m),
sd_dur_m = sd(total_m),
se_dur_m = std.error(total_m),
lci_dur_m = mean_dur_m - 1.96 * se_dur_m,
uci_dur_m = mean_dur_m + 1.96 * se_dur_m) %>%
ungroup()
}
### demographics by pair
trial_info_pair = df_trial_info %>%
group_by(pair) %>%
summarize_demographics() %>%
left_join(., df_condition_info) %>%
dplyr::select(pair, condition, total_s:trial_n)
### calculate mean experiment duration
mean_dur_cond = trial_info_pair %>%
group_by(condition) %>%
summarize_dur()
### summary statistics
trial_info_cond = df_trial_info %>%
group_by(condition) %>%
summarize_demographics() %>%
left_join(., mean_dur_cond) %>%
dplyr::select(condition, total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_cond
trial_info_pair_round = df_trial_info %>%
group_by(pair, round, round_n) %>%
summarize_demographics() %>%
left_join(., df_condition_info)
### calculate mean round duration
mean_dur_round = trial_info_pair_round %>%
group_by(round, round_n) %>%
summarize_dur()
trial_info_round = df_trial_info %>%
group_by(round, round_n) %>%
summarize_demographics() %>%
left_join(., mean_dur_round) %>%
dplyr::select(round, round_n, total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_round
### calculate mean duration by condition x round
mean_dur_cond_round = trial_info_pair_round %>%
group_by(condition, round, round_n) %>%
summarize_dur()
trial_info_cond_round = df_trial_info %>%
group_by(condition, round, round_n) %>%
summarize_demographics() %>%
left_join(., mean_dur_cond_round) %>%
dplyr::select(condition, round, round_n,
total_m, mean_dur_m:uci_dur_m,
everything(),
-ends_with("_s"))
trial_info_cond_round
### visibility condition: difference coding
h_cond = hypr(AO_Asym = AsymAV ~ AO,
Asym_Sym = SymAV ~ AsymAV,
levels = levels(df_trial_info$condition))
h_cond
## hypr object containing 2 null hypotheses:
## H0.AO_Asym: 0 = AsymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
##
## Call:
## hypr(AO_Asym = ~AsymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV",
## "AsymAV", "AO"))
##
## Hypothesis matrix (transposed):
## AO_Asym Asym_Sym
## SymAV 0 1
## AsymAV 1 -1
## AO -1 0
##
## Contrast matrix:
## AO_Asym Asym_Sym
## SymAV 1/3 2/3
## AsymAV 1/3 -1/3
## AO -2/3 -1/3
contrasts(df_trial_info$condition) = contr.hypothesis(h_cond)
### visibility condition: sum coding
h_cond = hypr(AO_Sym = SymAV ~ AO,
Asym_Sym = SymAV ~ AsymAV,
levels = levels(df_trial_info$condition))
h_cond
## hypr object containing 2 null hypotheses:
## H0.AO_Sym: 0 = SymAV - AO
## H0.Asym_Sym: 0 = SymAV - AsymAV
##
## Call:
## hypr(AO_Sym = ~SymAV - AO, Asym_Sym = ~SymAV - AsymAV, levels = c("SymAV",
## "AsymAV", "AO"))
##
## Hypothesis matrix (transposed):
## AO_Sym Asym_Sym
## SymAV 1 1
## AsymAV 0 -1
## AO -1 0
##
## Contrast matrix:
## AO_Sym Asym_Sym
## SymAV 1/3 1/3
## AsymAV 1/3 -2/3
## AO -2/3 1/3
contrasts(df_trial_info$condition_sum) = contr.hypothesis(h_cond)
### round
bacward_diff = hypr(R12 = R2 ~ R1,
R23 = R3 ~ R2,
R34 = R4 ~ R3,
R45 = R5 ~ R4,
R56 = R6 ~ R5,
levels = levels(df_trial_info$round))
bacward_diff
## hypr object containing 5 null hypotheses:
## H0.R12: 0 = R2 - R1
## H0.R23: 0 = R3 - R2
## H0.R34: 0 = R4 - R3
## H0.R45: 0 = R5 - R4
## H0.R56: 0 = R6 - R5
##
## Call:
## hypr(R12 = ~R2 - R1, R23 = ~R3 - R2, R34 = ~R4 - R3, R45 = ~R5 -
## R4, R56 = ~R6 - R5, levels = c("R1", "R2", "R3", "R4", "R5",
## "R6"))
##
## Hypothesis matrix (transposed):
## R12 R23 R34 R45 R56
## R1 -1 0 0 0 0
## R2 1 -1 0 0 0
## R3 0 1 -1 0 0
## R4 0 0 1 -1 0
## R5 0 0 0 1 -1
## R6 0 0 0 0 1
##
## Contrast matrix:
## R12 R23 R34 R45 R56
## R1 -5/6 -2/3 -1/2 -1/3 -1/6
## R2 1/6 -2/3 -1/2 -1/3 -1/6
## R3 1/6 1/3 -1/2 -1/3 -1/6
## R4 1/6 1/3 1/2 -1/3 -1/6
## R5 1/6 1/3 1/2 2/3 -1/6
## R6 1/6 1/3 1/2 2/3 5/6
contrasts(df_trial_info$round) = contr.hypothesis(bacward_diff)
bp_iconic_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=num_iconic_total, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Total N of iconic gestures per pair") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_iconic_by_cond)
bp_mean_iconic_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_num_iconic, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_num_iconic),
size = 2, shape = 4) +
scale_y_continuous(limits = c(0, 4), breaks = seq(0, 4, 1)) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean N of iconic gestures per trial") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_iconic_by_cond)
bp_mean_iconic_by_round_cond = ggplot(data=df_trial_info,
aes(x=round, y=num_iconic_gestures, fill=condition)) +
geom_boxplot(outlier.shape = NA,
alpha = 0.7) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 15), breaks = seq(0, 15, 5)) +
labs(x = "Round",
y = "Mean N of iconic gestures per trial") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
bp_mean_iconic_by_round_cond
The figure shows that the decline in iconic gestures across rounds does not follow a linear pattern. This suggests that log-transformed round may improve the model fit.
We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, which represents the expected number of iconic gestures in the SymAV condition, we will refer to the number of iconic gestures reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. Therefore, the expected number of iconic gestures per trial is \(4413 / (19 * 96) = 2.42\).
As Poisson regression models use a log-link, we need to specify priors on the log scale. Taking the log of 2.42, we get 0.88. In order to allow data to inform the posterior, we will use a normal prior with a mean of 0.88 and a relatively wide standard deviation of 3.
For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 2.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 2. For the correlation between the random effects, we set the prior to be LKJ(2).
For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean) in the SymAV condition. As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.
priors_rint = c(
prior(normal(0.88, 2), class = Intercept),
prior(normal(0, 2), class = b),
prior(normal(0, 2), class = sd))
priors_rslope = c(
prior(normal(0.88, 2), class = Intercept),
prior(normal(0, 2), class = b),
prior(normal(0, 2), class = sd),
prior(lkj(2), class = cor))
First, we will model the number of iconic gestures per trial with fixed effects for condition. As we expect no gestures in many trials based on previous study (Akamine et al., 2024), we will use a zero-inflated Poisson regression model with a log link function.
mp_iconic_cond_prior = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = poisson(),
prior = priors_rint,
sample_prior = "only",
data = df_trial_info,
file = "models/mp_iconic_cond_prior")
pp_check(mp_iconic_cond_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are similar to the observed data. The model generates a far wider range of data than the observed data, allowing the model to be informed by the data.
Next, we will fit the model to the observed data.
mp_iconic_cond = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = zero_inflated_poisson(),
prior = priors_rint,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mp_iconic_cond")
summary(mp_iconic_cond)
## Family: zero_inflated_poisson
## Links: mu = log; zi = identity
## Formula: num_iconic_gestures ~ 1 + condition + (1 | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.55 0.07 0.44 0.69 1.00 10772 22264
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.21 0.05 0.14 0.33 1.00 12904 27375
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 0.91 0.10 0.71 1.11 1.00 8817 16940
## conditionAO_Asym 0.29 0.21 -0.12 0.70 1.00 8349 15439
## conditionAsym_Sym -0.07 0.20 -0.48 0.33 1.00 8526 15481
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## zi 0.41 0.01 0.39 0.43 1.00 56110 51525
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(mp_iconic_cond)
bayestestR::hdi(mp_iconic_cond, ci = 0.89)
pp_check_sym = pp_check_each_condition(mp_iconic_cond, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(mp_iconic_cond, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(mp_iconic_cond, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
The posterior predictive check shows that the model generates data that are similar to the observed data. This suggests that the model is a good fit to the data.
Here, we will compare which regression model is suitable for the data, poisson or negative binomial. We will follow the procedure introduced in Winter & Burkner (2021).
First thing first, we will fit the negative binomial regression model to the data. Then we compare the predictive power of the models using loo_compare() function.
mnb_iconic_cond = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = negbinomial(),
prior = priors_rint,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond")
mnb_iconic_cond_zi = brm(num_iconic_gestures ~ 1 + condition +
(1|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rint,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_zi")
### loo compare
if (!file.exists("models/loo_pois_iconic.rds")){
pois_loo = loo(mp_iconic_cond)
saveRDS(pois_loo, file = "models/loo_pois_iconic.rds")
}
if (!file.exists("models/loo_nb_iconic.rds")){
nb_loo = loo(mnb_iconic_cond)
saveRDS(nb_loo, file = "models/loo_nb_iconic.rds")
}
if (!file.exists("models/loo_zinb_iconic.rds")){
zinb_loo = loo(mnb_iconic_cond_zi)
saveRDS(zinb_loo, file = "models/loo_zinb_iconic.rds")
}
pois_loo = readRDS("models/loo_pois_iconic.rds")
nb_loo = readRDS("models/loo_nb_iconic.rds")
zinb_loo = readRDS("models/loo_zinb_iconic.rds")
loo_compare(pois_loo, nb_loo, zinb_loo)
## elpd_diff se_diff
## mnb_iconic_cond 0.0 0.0
## mnb_iconic_cond_zi -0.9 1.5
## mp_iconic_cond -947.7 102.0
The loo_compare() function shows that the negative binomial regression has a larger predictive power (elpd_diff) and smaller standard error (se_diff) than the zero-inflated negative binomial regression and the Poisson regression. Therefore, we will use the negative binomial regression model for further analyses.
Next, we will run the negative binomial regression model to predict the number of iconic gestures per trial by condition and round. We will use the priors_round for the model.
We saw earlier that the decline in # iconic gestures doesn’t follow a linear pattern: The decline is stepper in the initial rounds and then plateaus afterwards. To capture this trend better without making 5 levels for Round (which is the case for bd coding), we will check if the model performance is comparable for backward-difference coded round and centered log-transformed round.
mnb_iconic_cond_round = brm(num_iconic_gestures ~ 1 + condition * round +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round")
mnb_iconic_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round_c")
mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round_log")
### loo compare
if (!file.exists("models/loo_nb_iconic_cond_round.rds")){
nb_cond_round_loo = loo(mnb_iconic_cond_round)
saveRDS(nb_cond_round_loo, file = "models/loo_nb_iconic_cond_round.rds")
}
if (!file.exists("models/loo_nb_iconic_cond_round_c.rds")){
nb_cond_round_c_loo = loo(mnb_iconic_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_iconic_cond_round_c.rds")
}
if (!file.exists("models/loo_nb_iconic_cond_round_log.rds")){
nb_cond_round_log_loo = loo(mnb_iconic_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_iconic_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/loo_nb_iconic_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_iconic_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)
## elpd_diff se_diff
## mnb_iconic_cond_round_log 0.0 0.0
## mnb_iconic_cond_round -13.9 5.3
## mnb_iconic_cond_round_c -38.2 12.8
The leave-one-out (LOO) Effect shows that centered log-transformed round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered round. This is in line with our observation that the decrease in the number of gestures does not follow a linear pattern, which explains why the centered round is disfavored here. Therefore, we will use the centered log round.
mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope,
sample_prior = "only",
data = df_trial_info,
file = "models/mnb_iconic_cond_round_log_prior")
pp_check(mnb_iconic_cond_round_log, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are similar to the observed data.
mnb_iconic_cond_round_log = brm(num_iconic_gestures ~ 1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_round,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/mnb_iconic_cond_round_log")
summary(mnb_iconic_cond_round_log)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_iconic_gestures ~ 1 + condition * log_round_c + (1 + log_round_c | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.85 0.10 0.68 1.07 1.00 15511
## sd(log_round_c) 0.41 0.06 0.31 0.54 1.00 24199
## cor(Intercept,log_round_c) 0.64 0.11 0.38 0.82 1.00 30037
## Tail_ESS
## sd(Intercept) 26870
## sd(log_round_c) 37632
## cor(Intercept,log_round_c) 44753
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.30 0.07 0.20 0.45 1.00 18463 33962
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -0.05 0.15 -0.34 0.25 1.00
## conditionAO_Asym 0.37 0.31 -0.24 0.98 1.00
## conditionAsym_Sym -0.16 0.31 -0.77 0.45 1.00
## log_round_c -1.37 0.07 -1.51 -1.23 1.00
## conditionAO_Asym:log_round_c 0.00 0.17 -0.33 0.34 1.00
## conditionAsym_Sym:log_round_c -0.10 0.17 -0.43 0.23 1.00
## Bulk_ESS Tail_ESS
## Intercept 12579 22655
## conditionAO_Asym 11906 21005
## conditionAsym_Sym 11866 21349
## log_round_c 18656 32465
## conditionAO_Asym:log_round_c 19212 31521
## conditionAsym_Sym:log_round_c 19007 32931
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 2.06 0.13 1.83 2.32 1.00 81667 50845
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(mnb_iconic_cond_round_log)
# bayestestR::hdi(mnb_iconic_cond_round_log, ci = 0.89)
Based on the model estimates, the condition does not have a significant effect on the number of iconic gestures produced per trial. However, the round has a significant effect on the number of iconic gestures produced per trial: the number of iconic gestures significantly decreases as the round progresses.
# p = plot_posterior(mnb_iconic_cond_round_log, interaction=T, xlim=1)
# p
p_direction(mnb_iconic_cond_round_log)
The probability of direction demonstrated significant effects for the round only.
### condition
#Sample the parameters of interest:
posterior <- as_draws_df(mnb_iconic_cond_round_log)
#Plot the prior-posterior update plot for the estimated parameters:
pp_update_intercept = ggplot(posterior) +
geom_density(aes(prior_Intercept), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_Intercept), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Intercept') +
theme_classic()
pp_update_b1 = ggplot(posterior) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAO_Asym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('AO --Asym') +
theme_classic()
pp_update_b2 = ggplot(posterior) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_conditionAsym_Sym), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Asym--Sym') +
theme_classic()
### round
pp_update_round = ggplot(posterior) +
geom_density(aes(prior_b), fill="steelblue", color="black",alpha=0.6) +
geom_density(aes(b_log_round_c), fill="#FC4E07", color="black",alpha=0.6) +
xlab('Centered log(round)') +
theme_classic()
gridExtra::grid.arrange(pp_update_intercept, pp_update_b1, pp_update_b2, pp_update_round)
pp_check_sym = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(mnb_iconic_cond_round_log, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
pp_check_r1 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(mnb_iconic_cond_round_log, df_trial_info, "R6")
gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)
emmeans(mnb_iconic_cond_round_log, pairwise ~ condition)$contrasts
## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV -0.160 -0.769 0.454
## SymAV - AO 0.214 -0.402 0.830
## AsymAV - AO 0.375 -0.240 0.981
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95
# emmeans(mnb_iconic_cond_round_log, pairwise ~ condition, level = 0.89)$contrasts
plot(conditional_effects(mnb_iconic_cond_round_log), ask = FALSE)
bp_mean_iconic_rate_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_iconic_per_100words, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .2,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_iconic_per_100words),
shape = 21, size = 3, fill = "white") +
scale_y_continuous(limits = c(0, 11), breaks = seq(0, 10, 2)) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean rate of iconic gestures") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_iconic_rate_by_cond)
pd = position_dodge(width = .75)
ggplot(data=df_trial_info,
aes(x=round, y=n_iconic_per_100words, fill=condition)) +
geom_boxplot(outlier.shape = NA,
alpha = 0.7) +
geom_point(data = trial_info_cond_round,
aes(y = mean_iconic_per_100words,
group = condition),
position = pd,
shape = 21, size = 2, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 22), breaks = seq(0, 20, 5)) +
labs(x = "Round",
y = "Mean rate of iconic gestures") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
As the unit of the dependent variable is different from the previous model, we will set different priors for the rate of iconic gestures per 100 words. We will use weakly informative priors for the regression each parameter. To specify priors for the intercept, we will refer to the rate of iconic gestures per 100 words reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of iconic gestures reported was 4413, which was collected from 19 dyads, each performing 96 trials. The total number of words produced was 71695 (28152 content) words. Therefore, the expected rate of iconic gestures per 100 words is \(4413 / (71695 / 100) = 6.16\) (log(6.16) = 1.82).
For the fixed effects, we will set unbiased priors with a mean of 0 and a standard deviation of 0.5.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).
For the models including the round as fixed effects (whether backward-difference coded or centered + log-transformed), the intercept will represent the mean expected number of iconic gestures (ground mean). As the meaning of the intercept doesn’t change when adding the round variable, we use the same prior for the intercept.
Note that we do not expect the rate of iconic gestures to change across rounds (i.e., we expect the number of words and gestures to decrease at an approximately same pace across the rounds). Also, it is common to set the mean of slopes to be 0 to avoid favoring any directions. Therefore we will set the mean of the prior for slope to 0.
priors_rslope_rate = c(
prior(normal(1.82, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_rate_zinb = c(
prior(normal(1.82, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))
In the previous section, we analyzed the number of iconic gestures produced per trial. However, it is common to analyze the rate of iconic gestures per 100 words to account for the differences in the length of the trials and speech rate. Here, we will include the log of speech rate (number of words / 100) as an exposure variable and analyze the rate of iconic gestures per 100 words by condition. Note that the syntax for the exposure variable is different from the Poisson regression model; for negative binomial regression, the exposure variable is included with the rate() function.
nb_iconic_rate_cond_round = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * round +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_iconic_rate_cond_round")
nb_iconic_rate_cond_round_c = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * round_c +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_iconic_rate_cond_round_c")
nb_iconic_rate_cond_round_log = brm(num_iconic_gestures | rate(n_words_100) ~
1 + condition * log_round_c +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_rate,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_iconic_rate_cond_round_log")
### loo compare
if (!file.exists("models/loo_nb_iconic_rate_cond_round.rds")){
nb_cond_round_loo = loo(nb_iconic_rate_cond_round)
saveRDS(nb_cond_round_loo, file = "models/loo_nb_iconic_rate_cond_round.rds")
}
if (!file.exists("models/loo_nb_iconic_rate_cond_round_c.rds")){
nb_cond_round_c_loo = loo(nb_iconic_rate_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_iconic_rate_cond_round_c.rds")
}
if (!file.exists("models/loo_nb_iconic_rate_cond_round_log.rds")){
nb_cond_round_log_loo = loo(nb_iconic_rate_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_iconic_rate_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/loo_nb_iconic_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_iconic_rate_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)
## elpd_diff se_diff
## nb_iconic_rate_cond_round_c 0.0 0.0
## nb_iconic_rate_cond_round -9.8 5.3
## nb_iconic_rate_cond_round_log -13.0 4.2
The leave-one-out (LOO) Effect shows that centered round provides a larger predictive power (elpd_diff) and smaller standard error (se_diff) compared to the backward-difference coded round or centered log-transformed round. Therefore, we will use the centered round.
zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
data = df_trial_info,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_iconic_rate_cond_round_c")
### loo compare
if (!file.exists("models/loo_zinb_iconic_rate_cond_round_c.rds")){
zinb_cond_round_c_loo = loo(zinb_iconic_rate_cond_round_c)
saveRDS(zinb_cond_round_c_loo, file = "models/loo_zinb_iconic_rate_cond_round_c.rds")
}
nb_cond_round_c_loo = readRDS("models/loo_nb_iconic_rate_cond_round_c.rds")
zinb_cond_round_c_loo = readRDS("models/loo_zinb_iconic_rate_cond_round_c.rds")
loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)
## elpd_diff se_diff
## zinb_iconic_rate_cond_round_c 0.0 0.0
## nb_iconic_rate_cond_round_c -13.7 6.2
The leave-one-out (LOO) Effect shows that zero-inflation model has a higher predictive power. As such, we will use zero-inflated negative binomial regression model.
zinb_iconic_rate_cond_round_c_prior = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
sample_prior = "only",
data = df_trial_info,
file = "models/zinb_iconic_rate_cond_round_c_prior")
pp_check(zinb_iconic_rate_cond_round_c_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model generates data that are somewhat similar to the observed data.
zinb_iconic_rate_cond_round_c = brm(num_iconic_gestures ~ 1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_rate_zinb,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_iconic_rate_cond_round_c")
model = zinb_iconic_rate_cond_round_c
summary(model)
## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = identity
## Formula: num_iconic_gestures ~ 1 + condition * round_c + offset(log(n_words_100)) + (1 + round_c | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.66 0.07 0.53 0.82 1.00 18554
## sd(round_c) 0.11 0.02 0.08 0.15 1.00 26924
## cor(Intercept,round_c) 0.71 0.11 0.46 0.87 1.00 38762
## Tail_ESS
## sd(Intercept) 34315
## sd(round_c) 44094
## cor(Intercept,round_c) 50332
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.10 0.03 0.06 0.16 1.00 26905 44443
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 1.33 0.10 1.13 1.54 1.00 12630
## conditionAO_Asym 0.28 0.22 -0.15 0.70 1.00 14119
## conditionAsym_Sym -0.06 0.22 -0.49 0.36 1.00 13682
## round_c -0.13 0.02 -0.17 -0.09 1.00 21391
## conditionAO_Asym:round_c 0.02 0.05 -0.07 0.11 1.00 24222
## conditionAsym_Sym:round_c -0.03 0.04 -0.12 0.06 1.00 22801
## Tail_ESS
## Intercept 23347
## conditionAO_Asym 25910
## conditionAsym_Sym 24753
## round_c 36599
## conditionAO_Asym:round_c 40686
## conditionAsym_Sym:round_c 36775
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 41.96 16.08 21.47 83.70 1.00 100817 52037
## zi 0.04 0.01 0.02 0.06 1.00 99853 44866
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
# bayestestR::hdi(model, ci = 0.89)
The coefficients show that the condition does not have a significant effect on the rate of iconic gestures per 100 words. However, there was a significant decrease in the rate of iconic gestures per 100 words across the rounds. This means that the number of iconic gestures decreased more than the number of words did across the rounds. A formal hypothesis testing will be performed later using Bayes factor.
# p = plot_posterior(model, interaction = T)
# p
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 1, 1.5)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(1.82, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))
fname = paste0("models/zinb_iconic_rate_cond_round_c_", prior_size[i])
fit = brm(num_iconic_gestures ~
1 + condition * round_c +
offset(log(n_words_100)) +
(1+round_c|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
# BF for sym - asym
h = hypothesis(fit, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
# BF for sym - ao
h = hypothesis(fit, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
# BF for round
h = hypothesis(fit, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round = c(bfs_round, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
# BF for sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])
# BF for sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])
# BF for round
h = hypothesis(model, "round_c = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round[3:5] = c(bf, bfs_round[3:4])
### make a df for BFs
df_bf = data.frame(size = prior_size,
sd = prior_sd,
ao_asym = bfs_cond_ao_asym,
asym_sym = bfs_cond_asym_sym,
round = bfs_round) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_asym", "asym_sym", "round"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Effect = factor(Effect,
levels = c("ao_asym", "asym_sym", "round"),
labels = c("AO-AsymAV", "AsymAV-SymAV", "Round")),
Predictor = ifelse(Effect == "round", "Round", "Visibility"),
BF10_log10 = log10(BF10))
df_bf %>% arrange(Effect, sd)
#### Plot BFs ####
ggplot(df_bf, aes(x = prior, y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor), scales="free_y") +
theme_bw(base_size = 12)+
theme(legend.position = "top")+
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("prior")
model = zinb_iconic_rate_cond_round_c
pp_check_sym = pp_check_each_condition(model, df_trial_info, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_trial_info, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_trial_info, "AO")
gridExtra::grid.arrange(pp_check_sym, pp_check_asym, pp_check_ao, ncol = 2)
pp_check_r1 = pp_check_each_round(model, df_trial_info, "R1")
pp_check_r2 = pp_check_each_round(model, df_trial_info, "R2")
pp_check_r3 = pp_check_each_round(model, df_trial_info, "R3")
pp_check_r4 = pp_check_each_round(model, df_trial_info, "R4")
pp_check_r5 = pp_check_each_round(model, df_trial_info, "R5")
pp_check_r6 = pp_check_each_round(model, df_trial_info, "R6")
gridExtra::grid.arrange(pp_check_r1, pp_check_r2, pp_check_r3, pp_check_r4, pp_check_r5, pp_check_r6, ncol = 3)
emmeans(model, pairwise ~ condition)$contrasts
## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV -0.0633 -0.497 0.355
## SymAV - AO 0.2134 -0.242 0.678
## AsymAV - AO 0.2760 -0.140 0.708
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95
emmeans(model, pairwise ~ condition, level = 0.89)$contrasts
## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV -0.0633 -0.405 0.286
## SymAV - AO 0.2134 -0.163 0.584
## AsymAV - AO 0.2760 -0.072 0.616
##
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.89
plot(conditional_effects(model), ask = FALSE)
Some experts in the field of statistics and causal inference have advised that researchers should build a causal model and examine which factors should be included and excluded from regression models if their aim is to infer the causal effects (e.g., McElreath, 2020, Pearl, 2010, Cinelli, Forney, & Pearl, 2020). Following this advice, we will build a causal model to infer the direct effect of speaker visibility on gestural alignment.
We assume the following causal model:
### causal model for gestural alignment rate
dag_gest <- dagitty('dag {
visibility -> gest_align
visibility -> n_words
visibility -> n_gestures
n_words -> lex_align
n_words -> n_gestures
round -> n_words
round -> lex_align
round -> gest_align
round -> n_gestures
n_gestures -> gest_align
lex_align -> gest_align
}')
plot(dag_gest)
### check the minimam adjustment set
print("Direct effect of visibility on gestural alignment rate")
## [1] "Direct effect of visibility on gestural alignment rate"
adjustmentSets(dag_gest, exposure = "visibility", outcome = "gest_align", effect = "direct")
## { lex_align, n_gestures, round }
## { n_gestures, n_words, round }
print("Direct effect of lexical alignment rate on gestural alignment rate")
## [1] "Direct effect of lexical alignment rate on gestural alignment rate"
adjustmentSets(dag_gest, exposure = "lex_align", outcome = "gest_align", effect = "direct")
## { n_gestures, round, visibility }
## { n_words, round }
print("Direct effect of round on gestural alignment rate")
## [1] "Direct effect of round on gestural alignment rate"
adjustmentSets(dag_gest, exposure = "round", outcome = "gest_align", effect = "direct")
## { lex_align, n_gestures, visibility }
The minimum adjustment set for estimating the direct causal effect of speaker visibility on gestural alignment rate is {lex_align, n_words, round}. Note that because dagitty didn’t find any adjustment set for the direct effect of visibility on lexical alignment rate when we expected bidirectional causation between lexical and gestural alignment, we assumed a unidirecitonal causation from lexical alignment to gestural alignment only in this model. This will be reversed in the causal model for lexical alignment rate, such that we assume a unidirectional causation from gestural alignment to lexical alignment.
In addition, we are also interested in the direct effect of lexical alignment rate on gestural alignment rate. The minimum adjustment set for is {visibility, n_gestures, round}.
As the minimum adjustment sets for the direct effects of visibility, lexical alignment rate, and round on gestural alignment rate are identical (i.e., {visibility, lex_align, n_gestures, round}), we can estimate the direct effect of these variables on gestural alignment rate with the following model:
\[ gest\_align \: | \: \text{rate}(n\_iconic\_gesture) \sim \\ \text{visibility} * \text{round} + \text{n_lexical_alignment} + \\ (1+\text{round} | \text{participant}) + (1 | \text{item}) \]
Just to be thorough, we will first analyze the raw frequency and rate
of gestural alignment. If you just want to see the analysis of the
proportion of gestural alignment, see [-–
bp_mean_gest_alignment_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_num_gest_align, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .4,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_num_gest_align),
size = 2, shape = 4) +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean gestural alignment count") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_gest_alignment_by_cond)
We will set priors for the intercept based on the expected number of gestural alignment reported in Akamine et al. (2024). In the paper, the authors analyzed data from 19 dyads involving in the same task as the current study but in co-present interaction. The total number of gestural alignment reported was 1086, which was collected from 19 dyads, each performing 96 trials. This leads to the expected number of gestural alignment is 1086 / (19 * 96) = 0.6 (log(0.6) = -0.5). As such, we will set the prior for the intercept to be normal with a mean of -0.5 and a standard deviation of 0.5.
For the fixed effects, we will set unbiased priors with a mean of 0. We set different SDs for each effect because the scale of each predictor is different. For N. lexical alignment and N. iconic gestures, we set a prior of Normal(0, 0.2). This means that if the mean N. lex align or N. iconic gesture increase by 1, we expect the mean N. gestural alignment to change by -0.6 to 0.9, with the most likely size of change to be 0. As for the visibility condition and the round, we set a prior of Normal(0. 0.5).
We will also specify a prior for the shape parameter to prevent the model from returning divergent transitions. We will set the prior to be normal with a mean of 0 and a wide standard deviation of 50.
For the standard deviation of the random effects, we set the prior to be normal with mean 0 and standard deviation 0.5. For the correlation between the random effects, we set the prior to be LKJ(2).
### pair
priors_rint_gest_align = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd))
priors_rslope_gest_align = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_gest_align_zinb = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.2), class = b, coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))
zinb_gest_align_prior = brm(num_gestural_align ~
1 + condition + round + num_lex_align + num_iconic_gestures +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_zinb,
data = df_trial_info,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/zinb_gest_align_prior")
pp_check(zinb_gest_align_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 4000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
zinb_align_cond_round = brm(num_gestural_align ~
1 + condition + round + num_lex_align + num_iconic_gestures +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_zinb,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_align_cond_round")
model = zinb_align_cond_round
summary(model)
## Family: zero_inflated_negbinomial
## Links: mu = log; shape = identity; zi = identity
## Formula: num_gestural_align ~ 1 + condition + round + num_lex_align + num_iconic_gestures + (1 + round | pair) + (1 | target)
## Data: df_trial_info (Number of observations: 4316)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 1.04 0.12 0.82 1.30 1.00 13258
## sd(roundR12) 0.62 0.15 0.33 0.93 1.00 24974
## sd(roundR23) 0.29 0.12 0.05 0.52 1.00 14676
## sd(roundR34) 0.19 0.12 0.01 0.46 1.00 12288
## sd(roundR45) 0.22 0.13 0.01 0.49 1.00 16214
## sd(roundR56) 0.16 0.12 0.01 0.45 1.00 28127
## cor(Intercept,roundR12) 0.33 0.21 -0.12 0.69 1.00 38507
## cor(Intercept,roundR23) 0.47 0.25 -0.13 0.85 1.00 36116
## cor(roundR12,roundR23) 0.06 0.27 -0.47 0.59 1.00 43615
## cor(Intercept,roundR34) 0.18 0.31 -0.48 0.72 1.00 44648
## cor(roundR12,roundR34) 0.10 0.31 -0.52 0.66 1.00 47546
## cor(roundR23,roundR34) 0.16 0.31 -0.48 0.72 1.00 41117
## cor(Intercept,roundR45) 0.15 0.31 -0.48 0.69 1.00 50329
## cor(roundR12,roundR45) 0.16 0.30 -0.46 0.70 1.00 45069
## cor(roundR23,roundR45) 0.20 0.31 -0.46 0.74 1.00 31625
## cor(roundR34,roundR45) 0.11 0.32 -0.54 0.69 1.00 35440
## cor(Intercept,roundR56) -0.08 0.33 -0.67 0.57 1.00 68954
## cor(roundR12,roundR56) 0.04 0.32 -0.58 0.64 1.00 62152
## cor(roundR23,roundR56) 0.04 0.33 -0.60 0.65 1.00 60123
## cor(roundR34,roundR56) 0.05 0.33 -0.59 0.66 1.00 50904
## cor(roundR45,roundR56) 0.03 0.33 -0.60 0.65 1.00 53360
## Tail_ESS
## sd(Intercept) 27474
## sd(roundR12) 27105
## sd(roundR23) 11704
## sd(roundR34) 20523
## sd(roundR45) 19835
## sd(roundR56) 30531
## cor(Intercept,roundR12) 46914
## cor(Intercept,roundR23) 32737
## cor(roundR12,roundR23) 47771
## cor(Intercept,roundR34) 48873
## cor(roundR12,roundR34) 53571
## cor(roundR23,roundR34) 49423
## cor(Intercept,roundR45) 48204
## cor(roundR12,roundR45) 49810
## cor(roundR23,roundR45) 44304
## cor(roundR34,roundR45) 48265
## cor(Intercept,roundR56) 53547
## cor(roundR12,roundR56) 53118
## cor(roundR23,roundR56) 57059
## cor(roundR34,roundR56) 56039
## cor(roundR45,roundR56) 56078
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.22 0.06 0.13 0.35 1.00 21997 37924
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -1.93 0.17 -2.26 -1.58 1.00 7722 15921
## conditionAO_Asym 0.32 0.29 -0.25 0.88 1.00 13178 24373
## conditionAsym_Sym 0.25 0.27 -0.29 0.79 1.00 11227 21534
## roundR12 1.54 0.16 1.23 1.85 1.00 31455 40797
## roundR23 -0.10 0.11 -0.32 0.10 1.00 30196 40085
## roundR34 -0.25 0.11 -0.47 -0.06 1.00 35515 38965
## roundR45 -0.31 0.12 -0.56 -0.08 1.00 39319 44558
## roundR56 -0.18 0.12 -0.43 0.05 1.00 58580 50167
## num_lex_align 0.05 0.02 0.00 0.09 1.00 68008 56486
## num_iconic_gestures 0.20 0.01 0.18 0.23 1.00 41732 49580
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 3.65 0.53 2.77 4.83 1.00 49170 51339
## zi 0.01 0.01 0.00 0.03 1.00 72219 39837
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
The coefficients show that the number of lexical alignment is a significant predictor of the number of gestural alignment.
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.05, 0.1, 0.3, 0.5)
bfs_lex_align = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(0.6, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b", coef = "num_lex_align"),
prior(normal(0, 0.2), class = b, coef = "num_iconic_gestures"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi),
prior(normal(0, 50), class = shape))
fname = paste0("models/zinb_align_cond_round_", prior_size[i])
fit = brm(num_gestural_align ~
1 + condition + round + num_lex_align + num_iconic_gestures +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors,
data = df_trial_info,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for N. lex alignment
h = hypothesis(fit, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align = c(bfs_lex_align, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.2, prior_sd[3:4])
### BF for N. lex alignment
h = hypothesis(model, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:5])
### make a df for BFs
df_bf_lex = data.frame(size = prior_size,
sd = prior_sd,
lex_align = bfs_lex_align) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("lex_align"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Effect = "N. lex align",
Predictor = "N. lex align")
#### Plot BFs ####
ggplot(filter(df_bf_lex),
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor)) +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
We can analyze the proportion of gestural alignment in two ways: (1) modeling the rate of gestural alignment per iconic gesture using a negative binomial regression model and (2) modeling the probability of gestural alignment using a zero-one-inflated beta regression model.
Key differences in the two models are that the negative binomial regression model assumes that the rate of gestural alignment is a count variable, while the zero-one-inflated beta regression model assumes that the proportion of gestural alignment is a continuous variable bounded between 0 and 1. Also, while negative binomial regression models the rate of events, the zero-one-inflated beta regression models the probability of events. In the case of the proportion of gestural alignment, two models should yield similar results, but it is important to note that the two models are conceptually different.
bp_mean_gest_alignment_prop_by_cond = ggplot(data=trial_info_pair,
aes(x=condition, y=mean_gest_align_prop, fill=condition)) +
geom_jitter(aes(color = pair),
size = 1, alpha = 1,
width = 0.07, height = 0) +
geom_boxplot(width = .3,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond,
aes(y = mean_gest_align_prop),
shape = 21, size = 3, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
labs(x = "Visibility",
y = "Mean gest align rate") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines"))
ggplotly(bp_mean_gest_alignment_prop_by_cond)
pd = position_dodge(width = .75)
bp_mean_gest_alignment_prop_by_round_cond =
ggplot(data=trial_info_pair_round,
aes(x=round_n, y=mean_gest_align_prop, fill = condition)) +
geom_jitter(aes(color = pair),
size = 0.5, alpha = 0.7,
width = 0.07, height = 0) +
geom_boxplot(width = .5,
outlier.shape = NA, alpha = 0.7) +
geom_point(data = trial_info_cond_round,
aes(y = mean_gest_align_prop),
position = pd,
shape = 21, size = 2, fill = "white") +
scale_fill_manual(values = c("#ED6B06", "#00786A", "darkgrey")) +
scale_y_continuous(limits = c(0, 1)) +
labs(x = "Round",
y = "Mean gest align rate") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "none",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
facet_grid(cols = vars(condition))
ggplotly(bp_mean_gest_alignment_prop_by_round_cond)
To model the proportion of gestural alignment, we need to remove trials where the number of iconic gestures is 0. This is because dividing any numbers by 0 results in an undefined value (NA) and prevents model from running.
df_gest_align_posreg_prop = df_trial_info %>%
filter(num_iconic_gestures > 0)
print(paste0("Number of rows before removing trials with no iconic gestures: ", nrow(df_trial_info)))
## [1] "Number of rows before removing trials with no iconic gestures: 4316"
print(paste0("Number of rows before after trials with no iconic gestures: ", nrow(df_gest_align_posreg_prop)))
## [1] "Number of rows before after trials with no iconic gestures: 2199"
print(paste0("Number of removed trials: ", nrow(df_trial_info) - nrow(df_gest_align_posreg_prop)))
## [1] "Number of removed trials: 2117"
We will set priors based on Akamine et al. (2024). As mentioned in the previous analysis, they detected 4413 iconic gestures and 1086 instances of gestural alignment. Dividing the number of gestural alignment by the number of iconic gestures gives 0.25 (1086/4413) (-1.39 on the log scale). This means that we expect 1 gestural alignment per 4 iconic gestures.
For the slope parameters, we set the mean of 0 with a SD of 0.5. This means that for example, we expect the mean difference between the AO and AsymAV conditions to range from -0.16 to 0.43.
### pair
priors_rint_gest_align_prop = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd))
priors_rslope_gest_align_prop = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 50), class = shape),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor))
priors_rslope_gest_align_prop_zinb = c(
prior(normal(-1.39, 0.5), class = Intercept),
prior(normal(0, 0.5), class = b),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 0.5), class = zi), # on the logit scale
prior(normal(0, 50), class = shape))
nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round")
nb_align_rate_cond_round_c = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round_c + num_lex_align +
(1+round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round_c")
nb_align_rate_cond_round_log = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * log_round_c + num_lex_align +
(1+log_round_c|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round_log")
### loo compare
if (!file.exists("models/loo_nb_align_rate_cond_round.rds")){
nb_cond_round_loo = loo(nb_align_rate_cond_round)
saveRDS(nb_cond_round_loo, file = "models/loo_nb_align_rate_cond_round.rds")
}
if (!file.exists("models/loo_nb_align_rate_cond_round_c.rds")){
nb_cond_round_c_loo = loo(nb_align_rate_cond_round_c)
saveRDS(nb_cond_round_c_loo, file = "models/loo_nb_align_rate_cond_round_c.rds")
}
if (!file.exists("models/loo_nb_align_rate_cond_round_log.rds")){
nb_cond_round_log_loo = loo(nb_align_rate_cond_round_log)
saveRDS(nb_cond_round_log_loo, file = "models/loo_nb_align_rate_cond_round_log.rds")
}
nb_cond_round_loo = readRDS("models/loo_nb_align_rate_cond_round.rds")
nb_cond_round_c_loo = readRDS("models/loo_nb_align_rate_cond_round_c.rds")
nb_cond_round_log_loo = readRDS("models/loo_nb_align_rate_cond_round_log.rds")
loo_compare(nb_cond_round_loo, nb_cond_round_c_loo, nb_cond_round_log_loo)
## elpd_diff se_diff
## nb_align_rate_cond_round 0.0 0.0
## nb_align_rate_cond_round_log -131.4 16.3
## nb_align_rate_cond_round_c -238.7 22.7
The leave-one-out (LOO) Effect shows that the backward-difference coded round provides a far larger predictive power (elpd_diff) and a far smaller standard error (se_diff) compared to the centered round or centered log-transformed round. Thus, we will use the bd coded round as a predictor for further analyses.
zinb_align_rate_cond_round = brm(num_gestural_align ~
1 + condition * round + num_lex_align +
offset(log(num_iconic_gestures)) +
(1+round|pair) + (1|target),
family = zero_inflated_negbinomial(),
prior = priors_rslope_gest_align_prop_zinb,
data = df_gest_align_posreg_prop,
sample_prior = T,
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/zinb_align_rate_cond_round")
### loo compare
if (!file.exists("models/loo_zinb_align_rate_cond_round.rds")){
zinb_cond_round_c_loo = loo(zinb_align_rate_cond_round)
saveRDS(zinb_cond_round_c_loo, file = "models/loo_zinb_align_rate_cond_round.rds")
}
nb_cond_round_c_loo = readRDS("models/loo_nb_align_rate_cond_round.rds")
zinb_cond_round_c_loo = readRDS("models/loo_zinb_align_rate_cond_round.rds")
loo_compare(nb_cond_round_c_loo, zinb_cond_round_c_loo)
## elpd_diff se_diff
## nb_align_rate_cond_round 0.0 0.0
## zinb_align_rate_cond_round -1.5 0.6
The loo comparsion shows that non-inflation model has a higher predictive power than the zero-inflation model. As such, we will use NB regression models for further analyses.
nb_gest_align_prop_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition + round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/nb_gest_align_prop_prior")
pp_check(nb_gest_align_prop_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
nb_align_rate_cond_round = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round")
model = nb_align_rate_cond_round
summary(model)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition * round + num_lex_align + (1 + round | pair) + (1 | target)
## Data: df_gest_align_posreg_prop (Number of observations: 2199)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.67 0.09 0.51 0.87 1.00 20129
## sd(roundR12) 0.64 0.12 0.43 0.90 1.00 32271
## sd(roundR23) 0.13 0.08 0.01 0.30 1.00 20697
## sd(roundR34) 0.08 0.07 0.00 0.25 1.00 31671
## sd(roundR45) 0.13 0.09 0.01 0.34 1.00 25468
## sd(roundR56) 0.13 0.10 0.00 0.36 1.00 35078
## cor(Intercept,roundR12) 0.03 0.21 -0.38 0.43 1.00 42101
## cor(Intercept,roundR23) 0.18 0.32 -0.49 0.74 1.00 63141
## cor(roundR12,roundR23) -0.22 0.30 -0.73 0.44 1.00 56062
## cor(Intercept,roundR34) -0.00 0.33 -0.63 0.63 1.00 86753
## cor(roundR12,roundR34) -0.09 0.32 -0.67 0.56 1.00 75896
## cor(roundR23,roundR34) 0.01 0.33 -0.63 0.63 1.00 65940
## cor(Intercept,roundR45) 0.04 0.33 -0.60 0.65 1.00 79946
## cor(roundR12,roundR45) 0.09 0.31 -0.53 0.65 1.00 71978
## cor(roundR23,roundR45) 0.01 0.33 -0.62 0.64 1.00 57760
## cor(roundR34,roundR45) -0.04 0.34 -0.66 0.61 1.00 53393
## cor(Intercept,roundR56) -0.06 0.34 -0.68 0.60 1.00 85023
## cor(roundR12,roundR56) 0.05 0.32 -0.57 0.64 1.00 79246
## cor(roundR23,roundR56) 0.00 0.33 -0.63 0.63 1.00 65196
## cor(roundR34,roundR56) -0.01 0.33 -0.64 0.63 1.00 58872
## cor(roundR45,roundR56) -0.02 0.33 -0.65 0.62 1.00 56595
## Tail_ESS
## sd(Intercept) 38159
## sd(roundR12) 47486
## sd(roundR23) 26398
## sd(roundR34) 34620
## sd(roundR45) 31341
## sd(roundR56) 32642
## cor(Intercept,roundR12) 50684
## cor(Intercept,roundR23) 50618
## cor(roundR12,roundR23) 49539
## cor(Intercept,roundR34) 52609
## cor(roundR12,roundR34) 53428
## cor(roundR23,roundR34) 57498
## cor(Intercept,roundR45) 53014
## cor(roundR12,roundR45) 54864
## cor(roundR23,roundR45) 57810
## cor(roundR34,roundR45) 56422
## cor(Intercept,roundR56) 53281
## cor(roundR12,roundR56) 54104
## cor(roundR23,roundR56) 56138
## cor(roundR34,roundR56) 58815
## cor(roundR45,roundR56) 58976
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.04 0.00 0.14 1.00 22309 20923
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -1.55 0.11 -1.77 -1.33 1.00 13026
## conditionAO_Asym 0.38 0.23 -0.08 0.84 1.00 16239
## conditionAsym_Sym 0.39 0.22 -0.06 0.83 1.00 12812
## roundR12 1.61 0.14 1.33 1.90 1.00 38765
## roundR23 0.24 0.08 0.09 0.39 1.00 54571
## roundR34 0.02 0.08 -0.14 0.19 1.00 58090
## roundR45 -0.12 0.10 -0.32 0.08 1.00 54198
## roundR56 -0.04 0.12 -0.27 0.19 1.00 64843
## num_lex_align -0.01 0.01 -0.04 0.02 1.00 100913
## conditionAO_Asym:roundR12 -0.07 0.28 -0.62 0.49 1.00 43120
## conditionAsym_Sym:roundR12 0.13 0.26 -0.39 0.64 1.00 36961
## conditionAO_Asym:roundR23 0.01 0.17 -0.32 0.34 1.00 57584
## conditionAsym_Sym:roundR23 -0.01 0.15 -0.30 0.29 1.00 55291
## conditionAO_Asym:roundR34 0.04 0.19 -0.33 0.40 1.00 55087
## conditionAsym_Sym:roundR34 0.15 0.16 -0.17 0.47 1.00 54834
## conditionAO_Asym:roundR45 0.24 0.22 -0.18 0.67 1.00 59488
## conditionAsym_Sym:roundR45 -0.14 0.19 -0.51 0.23 1.00 58188
## conditionAO_Asym:roundR56 0.19 0.26 -0.32 0.70 1.00 70505
## conditionAsym_Sym:roundR56 0.11 0.21 -0.30 0.51 1.00 72689
## Tail_ESS
## Intercept 26179
## conditionAO_Asym 30017
## conditionAsym_Sym 24995
## roundR12 48607
## roundR23 50032
## roundR34 50833
## roundR45 51064
## roundR56 54866
## num_lex_align 55847
## conditionAO_Asym:roundR12 52077
## conditionAsym_Sym:roundR12 47901
## conditionAO_Asym:roundR23 55257
## conditionAsym_Sym:roundR23 51500
## conditionAO_Asym:roundR34 55384
## conditionAsym_Sym:roundR34 54773
## conditionAO_Asym:roundR45 53403
## conditionAsym_Sym:roundR45 54752
## conditionAO_Asym:roundR56 55218
## conditionAsym_Sym:roundR56 57241
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 77.70 29.04 31.65 143.21 1.00 96781 51091
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
# bayestestR::hdi(model, ci = 0.89)
The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_asym = c()
bfs_cond_asym_sym = c()
bfs_round12 = c()
bfs_round23 = c()
bfs_round34 = c()
bfs_round45 = c()
bfs_round56 = c()
bfs_lex_align = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.39, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 50), class = shape))
fname = paste0("models/nb_align_rate_cond_round_", prior_size[i])
fit = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for visibility conditions
# sym - asym
h = hypothesis(fit, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym = c(bfs_cond_ao_asym, bf)
# sym - ao
h = hypothesis(fit, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym = c(bfs_cond_asym_sym, bf)
### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12 = c(bfs_round12, bf)
# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23 = c(bfs_round23, bf)
# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34 = c(bfs_round34, bf)
# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45 = c(bfs_round45, bf)
# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56 = c(bfs_round56, bf)
### BF for N. lex alignment
h = hypothesis(fit, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align = c(bfs_lex_align, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
### BF for visibility
# sym - asym
h = hypothesis(model, "conditionAO_Asym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_asym[3:5] = c(bf, bfs_cond_ao_asym[3:4])
# sym - ao
h = hypothesis(model, "conditionAsym_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_asym_sym[3:5] = c(bf, bfs_cond_asym_sym[3:4])
### BF for rounds
# R1 - R2
h = hypothesis(model, "roundR12 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round12[3:5] = c(bf, bfs_round12[3:4])
# R2 - R3
h = hypothesis(model, "roundR23 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round23[3:5] = c(bf, bfs_round23[3:4])
# R3 - R4
h = hypothesis(model, "roundR34 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round34[3:5] = c(bf, bfs_round34[3:4])
# R4 - R5
h = hypothesis(model, "roundR45 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round45[3:5] = c(bf, bfs_round45[3:4])
# R5 - R6
h = hypothesis(model, "roundR56 = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_round56[3:5] = c(bf, bfs_round56[3:4])
### BF for N. lex alignment
h = hypothesis(model, "num_lex_align = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_lex_align[3:5] = c(bf, bfs_lex_align[3:5])
### make a df for BFs
df_bf = data.frame(size = prior_size,
sd = prior_sd,
ao_asym = bfs_cond_ao_asym,
asym_sym = bfs_cond_asym_sym,
round12 = bfs_round12,
round23 = bfs_round23,
round34 = bfs_round34,
round45 = bfs_round45,
round56 = bfs_round56,
lex_align = bfs_lex_align) %>%
mutate(prior = paste0("N(0, ", sd, ")")) %>%
pivot_longer(cols = c("ao_asym", "asym_sym",
"round12", "round23", "round34", "round45", "round56",
"lex_align"),
names_to = "Effect",
values_to = "BF10") %>%
mutate(Predictor = ifelse(grepl("round", Effect), "Round",
ifelse(Effect == "lex_align", "N. lex align",
"Visibility")))
df_bf$Effect = recode(df_bf$Effect,
ao_asym = "AO--AsymAV",
asym_sym = "AsymAV--SymAV",
round12 = "R1--R2",
round23 = "R2--R3",
round34 = "R3--R4",
round45 = "R4--R5",
round56 = "R5--R6",
lex_align = "N. lex align")
#### Plot BFs ####
ggplot(filter(df_bf, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor)) +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
p_direction(model)
emmeans(model, pairwise ~ condition)$contrasts
## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.388 -0.0568 0.829
## SymAV - AO 0.770 0.2787 1.265
## AsymAV - AO 0.382 -0.0787 0.832
##
## Results are averaged over the levels of: round
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.95
emmeans(model, pairwise ~ condition, level = 0.89)$contrasts
## contrast estimate lower.HPD upper.HPD
## SymAV - AsymAV 0.388 0.028 0.744
## SymAV - AO 0.770 0.382 1.182
## AsymAV - AO 0.382 0.012 0.750
##
## Results are averaged over the levels of: round
## Point estimate displayed: median
## Results are given on the log (not the response) scale.
## HPD interval probability: 0.89
post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")
pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")
gridExtra::grid.arrange(pp_check_overall, pp_check_sym,
pp_check_asym, pp_check_ao,
ncol = 2)
Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.
plot(conditional_effects(model), ask=FALSE)
Here, we will run another model to compare the rate of gestural alignment between the SymAV and AO conditions.
nb_gest_align_prop_sum_prior = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum + round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = "only",
control = list(adapt_delta = 0.9,
max_treedepth = 20),
file = "models/nb_gest_align_prop_sum_prior")
pp_check(nb_gest_align_prop_sum_prior, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
The prior predictive check shows that the model expects fewer amount of 2, 3, and 4. This suggests that the zero-inflation prior may be too large or the mean for the intercept prior is too low. We will check the prior-posterior update plot and posterior predictive check to see if the model generates data that are similar to the observed data. If not, we will consider modifying the priors.
nb_align_rate_cond_round_sum = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors_rslope_gest_align_prop,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = "models/nb_align_rate_cond_round_sum")
model = nb_align_rate_cond_round_sum
summary(model)
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: num_gestural_align | rate(num_iconic_gestures) ~ 1 + condition_sum * round + num_lex_align + (1 + round | pair) + (1 | target)
## Data: df_gest_align_posreg_prop (Number of observations: 2199)
## Draws: 4 chains, each with iter = 20000; warmup = 2000; thin = 1;
## total post-warmup draws = 72000
##
## Multilevel Hyperparameters:
## ~pair (Number of levels: 45)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept) 0.67 0.09 0.51 0.87 1.00 18575
## sd(roundR12) 0.64 0.12 0.43 0.90 1.00 31125
## sd(roundR23) 0.12 0.08 0.01 0.30 1.00 20919
## sd(roundR34) 0.08 0.07 0.00 0.24 1.00 32022
## sd(roundR45) 0.12 0.09 0.00 0.33 1.00 26026
## sd(roundR56) 0.13 0.10 0.00 0.37 1.00 34132
## cor(Intercept,roundR12) 0.04 0.21 -0.38 0.44 1.00 40280
## cor(Intercept,roundR23) 0.18 0.32 -0.49 0.73 1.00 65536
## cor(roundR12,roundR23) -0.22 0.31 -0.73 0.44 1.00 59711
## cor(Intercept,roundR34) -0.01 0.34 -0.64 0.63 1.00 93535
## cor(roundR12,roundR34) -0.09 0.32 -0.67 0.56 1.00 78483
## cor(roundR23,roundR34) 0.01 0.33 -0.62 0.64 1.00 67466
## cor(Intercept,roundR45) 0.04 0.33 -0.60 0.65 1.00 80878
## cor(roundR12,roundR45) 0.09 0.31 -0.54 0.66 1.00 78177
## cor(roundR23,roundR45) 0.01 0.33 -0.62 0.64 1.00 57583
## cor(roundR34,roundR45) -0.04 0.33 -0.66 0.61 1.00 54005
## cor(Intercept,roundR56) -0.06 0.34 -0.68 0.59 1.00 93756
## cor(roundR12,roundR56) 0.05 0.32 -0.57 0.64 1.00 85093
## cor(roundR23,roundR56) -0.00 0.33 -0.63 0.63 1.00 67602
## cor(roundR34,roundR56) -0.01 0.33 -0.63 0.62 1.00 60265
## cor(roundR45,roundR56) -0.02 0.33 -0.64 0.62 1.00 58058
## Tail_ESS
## sd(Intercept) 34700
## sd(roundR12) 48295
## sd(roundR23) 26610
## sd(roundR34) 31291
## sd(roundR45) 32316
## sd(roundR56) 33988
## cor(Intercept,roundR12) 50279
## cor(Intercept,roundR23) 51923
## cor(roundR12,roundR23) 48590
## cor(Intercept,roundR34) 54422
## cor(roundR12,roundR34) 54685
## cor(roundR23,roundR34) 57533
## cor(Intercept,roundR45) 51517
## cor(roundR12,roundR45) 55396
## cor(roundR23,roundR45) 57358
## cor(roundR34,roundR45) 57922
## cor(Intercept,roundR56) 54259
## cor(roundR12,roundR56) 55731
## cor(roundR23,roundR56) 56162
## cor(roundR34,roundR56) 59440
## cor(roundR45,roundR56) 59143
##
## ~target (Number of levels: 16)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.06 0.03 0.00 0.14 1.00 21296 20498
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept -1.54 0.11 -1.77 -1.33 1.00
## condition_sumAO_Sym 0.65 0.23 0.19 1.11 1.00
## condition_sumAsym_Sym 0.27 0.23 -0.18 0.72 1.00
## roundR12 1.62 0.14 1.33 1.90 1.00
## roundR23 0.24 0.08 0.09 0.39 1.00
## roundR34 0.03 0.08 -0.13 0.19 1.00
## roundR45 -0.12 0.10 -0.32 0.08 1.00
## roundR56 -0.03 0.12 -0.27 0.20 1.00
## num_lex_align -0.01 0.01 -0.04 0.02 1.00
## condition_sumAO_Sym:roundR12 0.03 0.28 -0.52 0.57 1.00
## condition_sumAsym_Sym:roundR12 0.13 0.26 -0.39 0.65 1.00
## condition_sumAO_Sym:roundR23 0.01 0.17 -0.32 0.34 1.00
## condition_sumAsym_Sym:roundR23 -0.01 0.15 -0.30 0.29 1.00
## condition_sumAO_Sym:roundR34 0.16 0.18 -0.20 0.52 1.00
## condition_sumAsym_Sym:roundR34 0.15 0.16 -0.17 0.47 1.00
## condition_sumAO_Sym:roundR45 0.13 0.22 -0.29 0.56 1.00
## condition_sumAsym_Sym:roundR45 -0.15 0.19 -0.52 0.22 1.00
## condition_sumAO_Sym:roundR56 0.25 0.26 -0.26 0.76 1.00
## condition_sumAsym_Sym:roundR56 0.08 0.21 -0.33 0.49 1.00
## Bulk_ESS Tail_ESS
## Intercept 13062 24757
## condition_sumAO_Sym 13128 25189
## condition_sumAsym_Sym 11201 23088
## roundR12 38899 48092
## roundR23 54168 50887
## roundR34 61182 52654
## roundR45 58074 53364
## roundR56 66159 55148
## num_lex_align 95708 56583
## condition_sumAO_Sym:roundR12 43962 53222
## condition_sumAsym_Sym:roundR12 36727 47437
## condition_sumAO_Sym:roundR23 58050 53368
## condition_sumAsym_Sym:roundR23 60204 55326
## condition_sumAO_Sym:roundR34 60163 56713
## condition_sumAsym_Sym:roundR34 60303 55921
## condition_sumAO_Sym:roundR45 60850 56172
## condition_sumAsym_Sym:roundR45 60868 55720
## condition_sumAO_Sym:roundR56 70417 57299
## condition_sumAsym_Sym:roundR56 74686 56413
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 77.56 29.03 31.62 143.29 1.00 107322 49549
##
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
bayestestR::hdi(model)
# bayestestR::hdi(model, ci = 0.89)
The coefficients show that the proportion of gestural alignment was significantly higher in the SymAV condition than in the AO condition. Also, the rate of gestural alignment significantly increased from R1–R2 and R2–R3 and stabilized afterwards.
main_model = nb_align_rate_cond_round
model1 = zinb_align_cond_round
plot_posterior(main_model, model1, model)
### varying priors for sensitivity analysis
prior_size = c("xs", "s", "l", "xl")
prior_sd = c(0.1, 0.3, 0.7, 1)
bfs_cond_ao_sym = c()
for (i in 1:length(prior_sd)){
priors = c(
prior(normal(-1.39, 0.5), class = Intercept),
set_prior(paste0("normal(0,", prior_sd[i], ")"), class = "b"),
prior(normal(0, 0.5), class = sd),
prior(lkj(2), class = cor),
prior(normal(0, 50), class = shape))
fname = paste0("models/nb_align_rate_cond_round_sum_", prior_size[i])
fit = brm(num_gestural_align | rate(num_iconic_gestures) ~
1 + condition_sum * round + num_lex_align +
(1+round|pair) + (1|target),
family = negbinomial(),
prior = priors,
data = df_gest_align_posreg_prop,
sample_prior = T,
save_pars = save_pars(all = TRUE),
warmup = nwu, iter = niter,
control = list(adapt_delta = 0.9,
max_treedepth = 15),
file = fname)
### BF for visibility conditions
# BF for sym - ao
h = hypothesis(fit, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym = c(bfs_cond_ao_sym, bf)
}
### add BF for the main/medium model
prior_size[3:5] = c("m", prior_size[3:4])
prior_sd[3:5] = c(0.5, prior_sd[3:4])
### BF for visibility conditions
# sym - ao
h = hypothesis(model, "condition_sumAO_Sym = 0")
bf = 1 / abs(h$hypothesis$Evid.Ratio)
bfs_cond_ao_sym[3:5] = c(bf, bfs_cond_ao_sym[3:4])
### make a df for BFs
df_bf_temp = data.frame(size = prior_size,
sd = prior_sd,
Effect = "AO--SymAV",
BF10 = bfs_cond_ao_sym) %>%
mutate(prior = paste0("N(0, ", sd, ")"),
Predictor = "Visibility")
df_bf_new = df_bf %>%
filter(Effect != "N. lex align") %>%
rbind(df_bf_temp) %>%
rbind(df_bf_lex) %>%
mutate(Effect = factor(Effect,
levels = c("AO--SymAV", "AO--AsymAV", "AsymAV--SymAV",
"R1--R2", "R2--R3", "R3--R4", "R4--R5", "R5--R6",
"N. lex align")),
Predictor = factor(Predictor,
levels = c("Visibility", "Round", "N. lex align")))
df_bf_new %>% arrange(Effect, sd)
#### Plot BFs ####
ggplot(filter(df_bf_new, Effect!="R1--R2"), #exclude R1--R2 because BF is too huge
aes(x = factor(sd), y = BF10, group = Effect)) +
geom_hline(yintercept = 1, linetype="dashed") +
geom_point(aes(color=Effect)) +
geom_line(aes(color=Effect)) +
facet_wrap(vars(Predictor), scales = "free_x") +
theme_clean(base_size = 15) +
theme(axis.text.x = element_text(colour = "black", size = 14),
axis.text.y = element_text(colour = "black", size = 14),
axis.title = element_text(size = 15, face = 'bold'),
axis.title.x = element_text(vjust = -2),
axis.title.y = element_text(vjust = 2),
legend.position = "top",
legend.title=element_text(size=14),
legend.text=element_text(size=13),
strip.text = element_text(size = 15, face = 'bold'),
plot.background = element_blank(),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "lines")) +
scale_y_log10("Bayes factor (BF10)",
breaks = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100),
labels = c(0.001, 0.03, 0.01, 0.1, 0.33, 1, 3, 10, 30, 100)) +
xlab("SD for the prior")
p_direction(model)
post_sample = as_draws_df(model)
pp_update_plot(post_sample, model_type="nb")
pp_check_overall = pp_check(model, ndraws = 100, type = "bars") +
coord_cartesian(xlim = c(0, 20),
ylim = c(0, 3000))
pp_check_sym = pp_check_each_condition(model, df_gest_align_posreg_prop, "SymAV")
pp_check_asym = pp_check_each_condition(model, df_gest_align_posreg_prop, "AsymAV")
pp_check_ao = pp_check_each_condition(model, df_gest_align_posreg_prop, "AO")
gridExtra::grid.arrange(pp_check_overall, pp_check_sym,
pp_check_asym, pp_check_ao,
ncol = 2)
Although the model prediction is not perfect, this model had a higher predictive power than the negative binomial model. As such, we will use this model.
plot(conditional_effects(model), ask=FALSE)
cor = pcor.test(df_trial_info$num_lex_align, df_trial_info$num_gestural_align, df_trial_info$log_round_c)
cor
sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22631)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] rstan_2.32.6 StanHeaders_2.32.7 ppcor_1.1 MASS_7.3-58.1
## [5] doParallel_1.0.17 iterators_1.0.14 foreach_1.5.2 emmeans_1.10.7
## [9] tidybayes_3.0.7 bayestestR_0.15.2 brms_2.22.0 Rcpp_1.0.12
## [13] plotly_4.10.4 rsvg_2.6.2 DiagrammeRsvg_0.1 DiagrammeR_1.0.11
## [17] dagitty_0.3-4 ggh4x_0.3.1 ggthemes_5.1.0 hypr_0.2.8
## [21] plotrix_3.8-4 lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [25] dplyr_1.1.4 purrr_1.0.2 readr_2.1.5 tidyr_1.3.1
## [29] tibble_3.2.1 ggplot2_3.5.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.1-0 estimability_1.5.1 QuickJSR_1.1.3
## [4] rstudioapi_0.17.1 farver_2.1.1 svUnit_1.0.6
## [7] bit64_4.0.5 mvtnorm_1.2-4 bridgesampling_1.1-2
## [10] codetools_0.2-18 cachem_1.0.8 knitr_1.49
## [13] bayesplot_1.11.1 jsonlite_1.8.8 ggdist_3.3.2
## [16] compiler_4.2.2 httr_1.4.7 backports_1.4.1
## [19] Matrix_1.5-1 fastmap_1.1.1 lazyeval_0.2.2
## [22] cli_3.6.2 visNetwork_2.1.2 htmltools_0.5.8.1
## [25] tools_4.2.2 coda_0.19-4.1 gtable_0.3.6
## [28] glue_1.7.0 reshape2_1.4.4 posterior_1.6.1
## [31] V8_6.0.6 jquerylib_0.1.4 vctrs_0.6.5
## [34] nlme_3.1-160 crosstalk_1.2.1 insight_1.1.0
## [37] tensorA_0.36.2.1 xfun_0.53 ps_1.7.6
## [40] timechange_0.3.0 lifecycle_1.0.4 scales_1.3.0
## [43] vroom_1.6.5 ragg_1.3.0 hms_1.1.3
## [46] Brobdingnag_1.2-9 inline_0.3.21 RColorBrewer_1.1-3
## [49] yaml_2.3.8 curl_6.2.1 gridExtra_2.3
## [52] loo_2.8.0 sass_0.4.9 stringi_1.8.3
## [55] checkmate_2.3.1 pkgbuild_1.4.6 boot_1.3-28
## [58] cmdstanr_0.8.1 rlang_1.1.3 pkgconfig_2.0.3
## [61] systemfonts_1.0.6 matrixStats_1.3.0 distributional_0.5.0
## [64] pracma_2.4.4 evaluate_1.0.3 lattice_0.20-45
## [67] rstantools_2.4.0 htmlwidgets_1.6.4 labeling_0.4.3
## [70] processx_3.8.4 bit_4.0.5 tidyselect_1.2.1
## [73] plyr_1.8.9 magrittr_2.0.3 R6_2.6.1
## [76] generics_0.1.3 pillar_1.10.1 withr_3.0.2
## [79] datawizard_1.0.1 abind_1.4-8 crayon_1.5.3
## [82] arrayhelpers_1.1-0 tzdb_0.4.0 rmarkdown_2.29
## [85] grid_4.2.2 data.table_1.15.4 digest_0.6.35
## [88] xtable_1.8-4 textshaping_0.3.7 stats4_4.2.2
## [91] RcppParallel_5.1.7 munsell_0.5.1 viridisLite_0.4.2
## [94] bslib_0.9.0